source("choices.R")


### Read Data Start ###{{{
## Load from the file
read_raw_data <- function(file){
  raw_data <- read_csv(file, 
      col_types = cols(Debriefing_End_Time = col_time(format = "%H:%M:%S"), 
          Debriefing_Start_Time = col_time(format = "%H:%M:%S"), 
          EndDate = col_datetime(format = "%Y-%m-%d %H:%M:%S"), 
          GroupQ_End_Time = col_time(format = "%H:%M:%S"), 
          GroupQ_Start_Time = col_time(format = "%H:%M:%S"), 
          IdeologyQ_End_TIme = col_time(format = "%H:%M:%S"), 
          IdeologyQ_Start_Time = col_time(format = "%H:%M:%S"), 
          KnowledgeQ_End_Time = col_time(format = "%H:%M:%S"), 
          KnowledgeQ_Start_Time = col_time(format = "%H:%M:%S"), 
          RecordedDate = col_datetime(format = "%Y-%m-%d %H:%M:%S"), 
          StartDate = col_datetime(format = "%Y-%m-%d %H:%M:%S")), 
      na = "") # Set na="" when you download data from Qualtrics 
  names(raw_data) <- gsub(" ", "_", names(raw_data)) # Edit colnames
  names(raw_data) <- gsub("\\(|\\)", "", names(raw_data))

  print(paste("Observation num:", as.character(nrow(raw_data))))
  
  raw_data %>%
    select(-c(Status, IPAddress, RecipientEmail, RecipientFirstName,
              RecipientLastName, ExternalReference, LocationLatitude, LocationLongitude,
              nid, geoPC, geoPref, geoCity, Progress, Finished, DistributionChannel),
           -ends_with("First_Click"), -ends_with("Last_Click"),
           -ends_with("Page_Submit"), -ends_with("Click_Count"),
           -c(Duration_in_seconds, StartDate, EndDate, RecordedDate, ResponseId,
           Q_TotalDuration, KnowledgeQ_End_Time,
           GroupQ_Start_Time, GroupQ_End_Time, IdeologyQ_Start_Time,
           IdeologyQ_End_TIme, Debriefing_End_Time, Q1.1, Q16.4, Q16.5, Q17.4, Q17.5,
           Q4.1, Q8.1, Q15.6, Q15.3, Q15.4,
           # Unused questions
           Q2.2, Q2.3, Q5.3, Q5.4, Q5.6, Q5.7, 
           Q6.7, Q6.9, Q9.3, Q9.4, Q9.6, Q9.7, Q10.7, Q10.9, 
           Q12.2_1, Q12.2_2, Q12.2_3, Q13.4, Q14.4, Q14.5
           )) -> raw_data

  raw_data %>%
    mutate(Duration = Debriefing_Start_Time - KnowledgeQ_Start_Time) %>%
    select(-c(Debriefing_Start_Time, KnowledgeQ_Start_Time)) %>%
    rename(Completed = completed) -> raw_data

  return(raw_data)
}

clean_data <- function(raw_data, time, min_time=500, max_time=1500){
  ## Omit insufficient responses
  if(time == "GroupQ"){
    data <- raw_data %>%
      filter(Duration > get("min_time"), 
             Duration < get("max_time"), Completed == 1)

  }else{
    data <- raw_data %>%
        filter(Duration_in_seconds > get("min_time"), 
               Duration_in_seconds < get("max_time"), Completed == 1)
  }


  print(paste("Observation num:", as.character(nrow(data)),
              "/ Original:", as.character(nrow(raw_data)) ) )

  return(data)
}

### Read Data Finish}}}


### Global Variables {{{

choice_demography_all <- list(
                            c("Q15.3", "Q15.3_collapse", "Marriage"),
                            c("Q15.4", "Q15.4_collapse", "Children"),
                            c("Q15.7", "Q15.7_collapse_new", "Education")
                          )

choice_party_LDP <- list(c("Party_LDP"))
choice_party_Indp <- list(c("Party_Indp"))
choice_party_Others <- list(c("Party_Others"))
choice_party_Left <- list(c("Party_Left"))


choices_environment_ordered <- list(c("Q6.3", "Good", "DV1"),
                                    c("Q10.3", "Good", "DV2"))

choices_environment_binary <- list(c("Q6.3", "Good_collapse", "DV1"),
                                   c("Q10.3", "Good_collapse", "DV2"))

choice_rights <- list( 
                      Religion =  list(ordered = list(c("Q5.5", "Near_Reversed", "DV1"), 
                                                      # flipped values for consistency
                                                      c("Q9.5", "Near_Reversed", "DV2")),
                                       binary = list(c("Q5.5", "Near_collapse_new_Reversed", "DV1"),
                                                     c("Q9.5", "Near_collapse_new_Reversed", "DV2")),
                                       choices_rename = choices_good_eng),

                      Environment = list(ordered = choices_environment_ordered,
                                         binary = choices_environment_binary,
                                         choices_rename = choices_good_eng),

                      Transparency = list(ordered = list(c("Q6.4", "Good", "DV1"),
                                                         c("Q10.4", "Good", "DV2")),
                                          binary = list(c("Q6.4", "Good_collapse", "DV1"),
                                                        c("Q10.4", "Good_collapse", "DV2")),
                                          choices_rename = choices_good_eng),

                      Privacy =      list(ordered = list(c("Q6.5", "Good", "DV1"),
                                                         c("Q10.5", "Good", "DV2")),
                                          binary = list(c("Q6.5", "Good_collapse", "DV1"),
                                                        c("Q10.5", "Good_collapse", "DV2")),
                                          choices_rename = choices_good_eng),

                      Family =      list(ordered = list(c("Q6.6", "Good", "DV1"),
                                                         c("Q10.6", "Good", "DV2")),
                                          binary = list(c("Q6.6", "Good_collapse", "DV1"),
                                                        c("Q10.6", "Good_collapse", "DV2")),
                                          choices_rename = choices_good_eng),

                      VotingRights =      list(ordered = list(c("Q6.8", "Good", "DV1"),
                                                              c("Q10.8", "Good", "DV2")),
                                          binary = list(c("Q6.8", "Good_collapse", "DV1"),
                                                        c("Q10.8", "Good_collapse", "DV2")),
                                          choices_rename = choices_good_eng),

                      EmergencyDecree = list(ordered = list(c("Q7.4_1", "Good", "DV1"),
                                                           c("Q11.4_1", "Good", "DV2")),
                                          binary = list(c("Q7.4_1", "Good_collapse", "DV1"),
                                                        c("Q11.4_1", "Good_collapse", "DV2")),
                                          choices_rename = choices_good_eng),

                      EmergencyOrder = list(ordered = list(c("Q7.4_2", "Good", "DV1"),
                                                           c("Q11.4_2", "Good", "DV2")),
                                          binary = list(c("Q7.4_2", "Good_collapse", "DV1"),
                                                        c("Q11.4_2", "Good_collapse", "DV2")),
                                          choices_rename = choices_good_eng),

                      EmergencyDissolution = list(ordered = list(c("Q7.4_3", "Good", "DV1"),
                                                                c("Q11.4_3", "Good", "DV2")),
                                          binary = list(c("Q7.4_3", "Good_collapse", "DV1"),
                                                        c("Q11.4_3", "Good_collapse", "DV2")),
                                          choices_rename = choices_good_eng),


                      A9DefenceAmy =      list(ordered = list(c("Q7.3_1", "Good", "DV1"),
                                                              c("Q11.3_1", "Good", "DV2")),
                                          binary = list(c("Q7.3_1", "Good_collapse", "DV1"),
                                                        c("Q11.3_1", "Good_collapse", "DV2")),
                                          choices_rename = choices_good_eng),

                      A9IntCoop =         list(ordered = list(c("Q7.3_2", "Good", "DV1"),
                                                              c("Q11.3_2", "Good", "DV2")),
                                          binary = list(c("Q7.3_1", "Good_collapse", "DV1"),
                                                        c("Q11.3_1", "Good_collapse", "DV2")),
                                          choices_rename = choices_good_eng),

                      A9TerrResor =       list(ordered = list(c("Q7.3_3", "Good", "DV1"),
                                                              c("Q11.3_3", "Good", "DV2")),
                                          binary = list(c("Q7.3_3", "Good_collapse", "DV1"),
                                                        c("Q11.3_3", "Good_collapse", "DV2")),
                                          choices_rename = choices_good_eng)
                      )

figure_titles <- list(
                      FreedomOfSpeech = "Freedom of speech",
                      Equality = "Equality Under Law: Handicap",
                      Religion = "Separation of Church and State",
                      AmendmentRule = "Amendment Rule",
                      Emperor = "Emperor as HOS",

                      Environment = "Environment",
                      Transparency = "Government transparency",
                      Privacy = "Privacy",
                      Family = "Family Duty",
                      FlagAndAnthem = "Flag and Anthem",
                      VotingRights = "Foreigners' voting rights",
                      FiscalBalance = "Fiscal Balance",

                      EmergencyDecree = "Emergency: Decree",
                      EmergencyOrder = "Emergency: Follow orders",
                      EmergencyDissolution = "Emergency: No HR dissolution",
                      A9DefenceAmy = "Article 9: Defence Army",
                      A9IntCoop = "Article 9: International cooperation",
                      A9TerrResor = "Article 9: Territory and resources"

                     )


save_wide <- list(width = 12, height = 8)
save_square <- list(width = 9, height = 9)


### Global Variables }}}


### Utilities {{{

'%ni%' <- function(x,y)!('%in%'(x,y))


add_se <- function(model,name_f,name_x="Intercept"){
  #grab the standard error of the coefficients
  se_vec <- summary(model)$coefficients[,2]
  if(name_x=="Intercept"){
    #the stabdard error of the intercept
    se_x <- se_vec[1]
    #get the level-specific standard errors
    se_f <- se_vec[grep(name_f,names(se_vec))]
    se_f <- se_f[-grep(":",names(se_f))]
    #get the covariance between the intercept and the level-specific parameters
    vcov_f <- vcov(model)[grep(name_f,rownames(vcov(model))),grep(name_x,colnames(vcov(model)))]
    vcov_f <- vcov_f[-grep(":",names(vcov_f))]
    #the estimated average value at each level
    coef_f <- coef(model)[1]+coef(model)[names(vcov_f)]
  }
  else{
    #similar code for the case of another variable than the intercept
    se_x <- se_vec[name_x]
    se_f <- se_vec[grep(name_f,names(se_vec))]
    se_f <- se_f[grep(":",names(se_f))]
    vcov_f <- vcov(model)[grep(name_f,rownames(vcov(model))),grep(name_x,colnames(vcov(model)))][,1]
    vcov_f <- vcov_f[grep(":",names(vcov_f))]
    coef_f <- coef(model)[name_x]+coef(model)[names(vcov_f)]
  }
  #compute the summed SE
  se_f <- sqrt(se_x**2+se_f**2+2*vcov_f)
  #create the output dataframe
  out <- data.frame(Coef=coef_f,SE=se_f)
  return(out)
}


dummynumeric <- function(data, col) {
    # Can take numeric
    for(c in col) {
        idx <- which(names(data)==c)
        v <- data[[idx]]
        #stopifnot(class(v)=="factor")
        m <- matrix(0, nrow=nrow(data), ncol=length(sort(unique(v))))
        m[cbind(seq_along(v), as.integer(v))]<-1
        colnames(m) <- paste(c, sort(unique(v)), sep="_")
        r <- data.frame(m)
        if ( idx>1 ) {
            r <- cbind(data[1:(idx-1)],r)
        }
        if ( idx<ncol(data) ) {
            r <- cbind(r, data[(idx+1):ncol(data)])
        }
        data <- r
    }
    data
}

###}}}



### Analysis ###{{{

func_pass <- function(...){
  # Placeholder
  message("Passed")
}

analysis_settings_h <- function(){
  print("analysis_setings = list(foldername, sample, qid, cov_type, interaction_type, method)")
  print("OR")
  print("analysis_settings = list(foldername, sample, qid, file_name)")
}


save_obj <- function(obj, analysis_settings){

  savename <- make_filename(analysis_settings, ".obj")
  saveRDS(obj, file=savename)

}

check_obj <- function(analysis_settings){
  
  filename <- make_filename(analysis_settings, ".obj")

  if(file.exists(paste0(filename))){
    message(paste0("File alredy exists, we skip analysis: ", filename))
    return(readRDS(filename)$returned)  
  }
  
  return(NULL)

}

analysis <- function(data, change_choices_set, 
                     analysis_func=func_pass, figure_func=func_pass, 
                     analysis_settings=list(), figure_settings=list(), 
                     save_settings=list(),
                     return_res = F)
{
  savefilename <- make_filename(analysis_settings)
  exist_figure <- check_figure_existence(savefilename)
  if(exist_figure)
    message(paste0("File alredy exists, we skip it: ", savefilename)) ; return

  tempdata <- data_change_choices(data, change_choices_set)

  if(is.function(analysis_func)){
    res <- check_obj(analysis_settings)

    if(is.null(res))
      res <- analysis_func(tempdata, analysis_settings)

    p <- figure_func(res, analysis_settings, figure_settings)
    p <- save_figure(p, analysis_settings, save_settings)  
  }else{
    stop("Check `analysis_func` (seems it's not a function)")  
  }

  if(return_res)
    return(res)


}


analysis_run_Stan <- function(model, stan_data, chains=5, iter=1000, seed=919){

  return(sampling(model, data=stan_data, chains=chains, 
                  iter=iter, seed=seed,
                  control = list(max_treedepth = 12)))

}


analysis_diagnosis_Stan <- function(resStan, analysis_settings){
  filename <- make_filename(analysis_settings)
  filename <- str_replace(filename, "\\.pdf", "_Rhat\\.pdf")

  p <- stan_rhat(resStan)
  save_figure(p, analysis_settings, save_settings=list(width=6, height=6), filename)

}


analysis_bayes_ordered <- function(tempdata, analysis_settings)
{
  # BayesSimpla: no Cov, no Focus

  tempdata %>% 
    mutate(depend_var = if_else(Group==0, DV1, DV2)) %>%
    filter(depend_var != 0) -> temp

  stan_data <- list(
                    N = nrow(temp),
                    K = length(unique(temp$depend_var)),
                    treatment = temp %>% pull(Group),
                    Y = temp %>% pull(depend_var)
                   )

  model <- ordered_logit_noCov

  resStan <- analysis_run_Stan(model, stan_data,
                               iter=analysis_settings$stan_iteration)

  save_obj(list(returned=resStan), analysis_settings)
  analysis_diagnosis_Stan(resStan, analysis_settings)

  return(resStan)
}


analysis_bayes_ordered_focus <- function(tempdata, analysis_settings)
{
  # No cov but focus

  focus_var_name <- analysis_settings$focus_var_name

  tempdata %>% 
    mutate(depend_var = if_else(Group==0, DV1, DV2)) %>%
    filter(depend_var != 0) -> temp

  stan_data <- list(
                    N = nrow(temp),
                    K = length(unique(temp$depend_var)),
                    F = length(unique(temp %>% pull(focus_var_name))),
                    focus = temp %>% pull(focus_var_name),
                    treatment = temp %>% pull(Group),
                    Y = temp %>% pull(depend_var)
                   )

  model <- ordered_logit2_noCov

  resStan <- analysis_run_Stan(model, stan_data,
                               iter=analysis_settings$stan_iteration)

  save_obj(list(returned=resStan), analysis_settings)
  analysis_diagnosis_Stan(resStan, analysis_settings)

  return(resStan)

}


analysis_bayes_cov <- function(tempdata, analysis_settings){
  # Cov but no focus

  tempdata %>% 
    mutate(depend_var = if_else(Group==0, DV1, DV2)) %>%
    filter(depend_var != 0) -> temp
  
  cov_names <- analysis_settings$cov_names
  covariates <- temp %>% select(cov_names) %>% as.matrix()

  stan_data <- list(
                    N = nrow(temp),
                    K = length(unique(temp$depend_var)),
                    M = ncol(covariates),
                    covariates = covariates,
                    treatment = temp %>% pull(Group),
                    Y = temp %>% pull(depend_var)
                   )

  model <- ordered_logit_noCov

  resStan <- analysis_run_Stan(model, stan_data,
                               iter=analysis_settings$stan_iteration)

  save_obj(list(returned=resStan), analysis_settings)
  analysis_diagnosis_Stan(resStan, analysis_settings)

  return(resStan)


}


analysis_bayes_cov_focus <- function(tempdata, analysis_settings){
  # Cov and focus

  focus_var_name <- analysis_settings$focus_var_name

  tempdata %>% 
    mutate(depend_var = if_else(Group==0, DV1, DV2)) %>%
    filter(depend_var != 0) -> temp

  cov_names <- analysis_settings$cov_names
  covariates <- temp %>% select(cov_names) %>% as.matrix()

  stan_data <- list(
                    N = nrow(temp),
                    K = length(unique(temp$depend_var)),
                    M = ncol(covariates),
                    F = length(unique(temp %>% pull(focus_var_name))),
                    covariates = covariates,
                    focus = temp %>% pull(focus_var_name),
                    treatment = temp %>% pull(Group),
                    Y = temp %>% pull(depend_var)
                   )

  model <- ordered_logit2

  resStan <- analysis_run_Stan(model, stan_data,
                               iter=analysis_settings$stan_iteration)

  save_obj(list(returned=resStan), analysis_settings)
  analysis_diagnosis_Stan(resStan, analysis_settings)

  return(resStan)


}



analysis_freq_parametric <- function(tempdata, analysis_settings)
{
  formulas <- analysis_settings$formula
  focus_var_name <- analysis_settings$focus_var_name
  interaction_name <- analysis_settings$interaction_name

  # Full ordered
  tempdata %>% 
    mutate(depend_var = if_else(Group==0, DV1, DV2)) %>%
    filter(depend_var != 0) %>%
    mutate(Outcome = 6 - depend_var) %>%
    rename(Treatment = Group) -> temp

  message("In Freq Ordered, Values are reversed (Good = 5, Bad = 1)")

  res_lm1 <- lm(as.formula(formulas$FreqOrdered), data=temp)

  if(is.null(focus_var_name))
    res_lm1_broom <- broom::tidy(res_lm1)

  if(!is.null(focus_var_name)){
    res_lm1_broom <- broom::tidy(res_lm1) %>% filter(term=="Treatment") %>%
                      select(term, estimate, std.error)
    interact <- as_tibble(add_se(res_lm1, name_f=focus_var_name, name_x="Treatment")) %>%
                  mutate(term = "Treatment") %>%
                  rename(estimate = Coef,
                         std.error = SE)

    res_lm1_broom <- bind_rows(res_lm1_broom, interact)
    res_lm1_broom$Group <- interaction_name
                           
  }


  # Collapsed (Difference-in-means)
  tempdata %>% 
    mutate(depend_var = if_else(Group==0, DV1, DV2)) %>%
    mutate(Outcome = case_when(
                                  depend_var == 1 ~ 1,
                                  depend_var == 2 ~ 1,
                                  depend_var == 3 ~ -1,
                                  depend_var == 4 ~ 0,
                                  depend_var == 5 ~ 0,
                                  depend_var == 0 ~ -1,
                               )) %>%
    filter(Outcome >= 0) %>%
    rename(Treatment = Group) -> temp

  res_lm2 <- lm(as.formula(formulas$FreqCollapsed), data=temp)

  if(is.null(focus_var_name))
    res_lm2_broom <- broom::tidy(res_lm2)

  if(!is.null(focus_var_name)){
    res_lm2_broom <- broom::tidy(res_lm2) %>% filter(term=="Treatment") %>%
                      select(term, estimate, std.error)
    interact <- as_tibble(add_se(res_lm2, name_f=focus_var_name, name_x="Treatment")) %>%
                  mutate(term = "Treatment") %>%
                  rename(estimate = Coef,
                         std.error = SE)

    res_lm2_broom <- bind_rows(res_lm2_broom, interact)
    res_lm2_broom$Group <- interaction_name
                           
  }


  res_lm1_broom$Model <- "5 choices"
  res_lm2_broom$Model <- "2 choices"

  res <- bind_rows(res_lm1_broom, res_lm2_broom)
    
  save_obj(list(raw1=res_lm1, raw2=res_lm2, returned=res),
           analysis_settings)

  return(res)

}



analysis_freq_nonparametric <- function(tempdata, analysis_settings){

  choices_name <- analysis_settings$choices_rename
  formula <- analysis_settings$formula$FreqNonparametric
  focus_var_name <- analysis_settings$focus_var_name
  interaction_name <- analysis_settings$interaction_name


  # 1
  tempdata %>% 
    mutate(depend_var = if_else(Group==0, DV1, DV2)) %>%
    filter(depend_var != 0) %>%
    mutate(depend_var = case_when(
                                  depend_var == 1 ~ 1,
                                  depend_var == 2 ~ 0,
                                  depend_var == 3 ~ 0,
                                  depend_var == 4 ~ 0,
                                  depend_var == 5 ~ 0
                                 )) %>%
    rename(Treatment = Group, Outcome=depend_var) -> temp

  res_lm1 <- lm(as.formula(formula), data=temp)

  if(is.null(focus_var_name))
    res_lm1_broom <- broom::tidy(res_lm1)

  if(!is.null(focus_var_name)){
    res_lm1_broom <- broom::tidy(res_lm1) %>% filter(term=="Treatment") %>%
                      select(term, estimate, std.error)
    interact <- as_tibble(add_se(res_lm1, name_f=focus_var_name, name_x="Treatment")) %>%
                  mutate(term = "Treatment") %>%
                  rename(estimate = Coef,
                         std.error = SE)

    res_lm1_broom <- bind_rows(res_lm1_broom, interact)
    res_lm1_broom$Group <- interaction_name
                           
  }

  res_lm1_broom$term <- gsub("Treatment", choices_name[1], res_lm1_broom$term)


  # 2
  tempdata %>% 
    mutate(depend_var = if_else(Group==0, DV1, DV2)) %>%
    filter(depend_var != 0) %>%
    mutate(depend_var = case_when(
                                  depend_var == 1 ~ 0,
                                  depend_var == 2 ~ 1,
                                  depend_var == 3 ~ 0,
                                  depend_var == 4 ~ 0,
                                  depend_var == 5 ~ 0
                                 )) %>%
    rename(Treatment = Group, Outcome=depend_var) -> temp

  res_lm2 <- lm(as.formula(formula), data=temp)

  if(is.null(focus_var_name))
    res_lm2_broom <- broom::tidy(res_lm2)

  if(!is.null(focus_var_name)){
    res_lm2_broom <- broom::tidy(res_lm2) %>% filter(term=="Treatment") %>%
                      select(term, estimate, std.error)
    interact <- as_tibble(add_se(res_lm2, name_f=focus_var_name, name_x="Treatment")) %>%
                  mutate(term = "Treatment") %>%
                  rename(estimate = Coef,
                         std.error = SE)

    res_lm2_broom <- bind_rows(res_lm2_broom, interact)
    res_lm2_broom$Group <- interaction_name
                           
  }


  res_lm2_broom$term <- gsub("Treatment", choices_name[2], res_lm2_broom$term)


  # 3
  tempdata %>% 
    mutate(depend_var = if_else(Group==0, DV1, DV2)) %>%
    filter(depend_var != 0) %>%
    mutate(depend_var = case_when(
                                  depend_var == 1 ~ 0,
                                  depend_var == 2 ~ 0,
                                  depend_var == 3 ~ 1,
                                  depend_var == 4 ~ 0,
                                  depend_var == 5 ~ 0
                                 )) %>%
    rename(Treatment = Group, Outcome=depend_var) -> temp

  res_lm3 <- lm(as.formula(formula), data=temp)

  if(is.null(focus_var_name))
    res_lm3_broom <- broom::tidy(res_lm3)

  if(!is.null(focus_var_name)){
    res_lm3_broom <- broom::tidy(res_lm3) %>% filter(term=="Treatment") %>%
                      select(term, estimate, std.error)
    interact <- as_tibble(add_se(res_lm3, name_f=focus_var_name, name_x="Treatment")) %>%
                  mutate(term = "Treatment") %>%
                  rename(estimate = Coef,
                         std.error = SE)

    res_lm3_broom <- bind_rows(res_lm3_broom, interact)
    res_lm3_broom$Group <- interaction_name
                           
  }


  res_lm3_broom$term <- gsub("Treatment", choices_name[3], res_lm3_broom$term)


  # 4
  tempdata %>% 
    mutate(depend_var = if_else(Group==0, DV1, DV2)) %>%
    filter(depend_var != 0) %>%
    mutate(depend_var = case_when(
                                  depend_var == 1 ~ 0,
                                  depend_var == 2 ~ 0,
                                  depend_var == 3 ~ 0,
                                  depend_var == 4 ~ 1,
                                  depend_var == 5 ~ 0
                                 )) %>%
    rename(Treatment = Group, Outcome=depend_var) -> temp

  res_lm4 <- lm(as.formula(formula), data=temp)

  if(is.null(focus_var_name))
    res_lm4_broom <- broom::tidy(res_lm4)

  if(!is.null(focus_var_name)){
    res_lm4_broom <- broom::tidy(res_lm4) %>% filter(term=="Treatment") %>%
                      select(term, estimate, std.error)
    interact <- as_tibble(add_se(res_lm4, name_f=focus_var_name, name_x="Treatment")) %>%
                  mutate(term = "Treatment") %>%
                  rename(estimate = Coef,
                         std.error = SE)

    res_lm4_broom <- bind_rows(res_lm4_broom, interact)
    res_lm4_broom$Group <- interaction_name
                           
  }


  res_lm4_broom$term <- gsub("Treatment", choices_name[4], res_lm4_broom$term)


  # 5
  tempdata %>% 
    mutate(depend_var = if_else(Group==0, DV1, DV2)) %>%
    filter(depend_var != 0) %>%
    mutate(depend_var = case_when(
                                  depend_var == 1 ~ 0,
                                  depend_var == 2 ~ 0,
                                  depend_var == 3 ~ 0,
                                  depend_var == 4 ~ 0,
                                  depend_var == 5 ~ 1
                                 )) %>%
    rename(Treatment = Group, Outcome=depend_var) -> temp

  res_lm5 <- lm(as.formula(formula), data=temp)

  if(is.null(focus_var_name))
    res_lm5_broom <- broom::tidy(res_lm5)

  if(!is.null(focus_var_name)){
    res_lm5_broom <- broom::tidy(res_lm5) %>% filter(term=="Treatment") %>%
                      select(term, estimate, std.error)
    interact <- as_tibble(add_se(res_lm5, name_f=focus_var_name, name_x="Treatment")) %>%
                  mutate(term = "Treatment") %>%
                  rename(estimate = Coef,
                         std.error = SE)

    res_lm5_broom <- bind_rows(res_lm5_broom, interact)
    res_lm5_broom$Group <- interaction_name
                           
  }


  res_lm5_broom$term <- gsub("Treatment", choices_name[5], res_lm5_broom$term)


  res <- bind_rows(res_lm1_broom, res_lm2_broom, res_lm3_broom,
                   res_lm4_broom, res_lm5_broom)


  save_obj(list(raw1=res_lm1,
                raw2=res_lm2,
                raw3=res_lm3,
                raw4=res_lm4,
                raw5=res_lm5,
                returned=res),
           analysis_settings)

  return(res)

}


## Run main analysis

main_analysis <- function(data,
                          choice_rights, 
                          cov_types, cov_types_names,
                          interaction_types, interaction_names,
                          analysis_functions, figure_functions
                         )

{

  for(right in names(choice_rights)){

  for(j in 1:length(cov_types)){
  cov_type <- cov_types[[j]]


  for(k in 1:length(interaction_types)){
  interaction_type <- interaction_types[[k]]


  # Starting each loop element
  analysis_todo <- c()

  if(names(cov_types)[j] != "None" & names(interaction_types)[k] == "None"){
    # cov, no itneraction -> BayesCov

    interaction_type <- interaction_types[[k]]
    # analysis_todo <- c(analysis_todo, "BayesCov", "FreqParametric", "FreqNonparametric")
    analysis_todo <- c(analysis_todo, "BayesCov")
    choices <- choice_rights[[right]]$ordered
    choices_rename <- choice_rights[[right]]$choices_rename
    focus_var_name <- NULL
    focus_names <- NULL
    interaction_name <- NULL

    cov_type <- cov_types[[j]]
    cov_names <- unlist(lapply(cov_type, function(x){x[3]}))

    cov_term <- paste0("Treatment + ", paste(cov_names, collapse=" + "))
    formulas_list <- list(
                            FreqOrdered = paste0("Outcome ~ ", cov_term),
                            FreqCollapsed = paste0("Outcome ~ ", cov_term),
                            FreqNonparametric = paste0("Outcome ~ ", cov_term)
                         )
  }


  if(names(cov_types)[j] != "None" & names(interaction_types)[k] != "None"){
    # cov and iteraction -> BayesCovFocus

    interaction_type <- interaction_types[[k]]
    # analysis_todo <- c(analysis_todo, "BayesCovFocus", "FreqParametric", "FreqNonparametric")
    analysis_todo <- c(analysis_todo, "BayesCovFocus")
    choices <- choice_rights[[right]]$ordered
    choices_rename <- choice_rights[[right]]$choices_rename

    focus_var_name <- interaction_type[[1]][3]
    focus_names <- interaction_type[[1]][3]
    interaction_name <- interaction_names[[k]][[1]]

    cov_type <- cov_types[[j]]
    cov_names <- unlist(lapply(cov_type, function(x){x[3]}))
    cov_term <- paste(cov_names, collapse=" + ")

    interaction_term <- paste0("Treatment*", interaction_type[[1]][3], " + ", cov_term)
    formulas_list <- list(
                            FreqOrdered = paste0("Outcome ~ ", interaction_term),
                            FreqCollapsed = paste0("Outcome ~ ", interaction_term),
                            FreqNonparametric = paste0("Outcome ~ ", interaction_term)
                         )
  }


  for(an in 1:length(analysis_todo)){
    analysis_name <- analysis_todo[an]

    analysis(
              data,
              change_choices_set = c(Group=list(c("Group", "Group_new", "Group")),
                                     DV1=choices[1], DV2=choices[2],
                                     cov_type, interaction_type),
              analysis_settings = list(paste0(right, "/"), sample_name, right, 
                                       names(cov_types)[j], names(interaction_types)[k], 
                                       analysis_name,
                                       stan_iteration=stan_iteration,
                                       choices_rename=choices_rename,
                                       formula=formulas_list,
                                       focus_var_name=focus_var_name,
                                       interaction_name = interaction_name,
                                       cov_names = cov_names
                                      ),
              analysis_func = analysis_functions[[analysis_name]],
              figure_func = figure_functions[[analysis_name]],
              figure_settings = list(title=figure_titles[[right]], choices_rename=choices_rename,
                                     focus_var_name=focus_var_name, focus_names=focus_names),
              save_settings = save_square,
              return_res = F
            )
  }

  }

  }


  }

}


### }}}


### Figure{{{

make_filename <- function(analysis_settings, extention=".pdf"){
  # Make a filename from the analysis settings
  # analysis_setings = list(foldername, sample, qid, cov_type, interaction_type, method)
  # OR
  # analysis_settings = list(foldername, sample, qid, file_name)
  #
  # return full_path

  if(!is.list(analysis_settings))
    stop("Provide list to make_filename()")


  foldername <- paste0(savefolder_base, analysis_settings[[1]])
  if(!file.exists(foldername)){
    message("Figure folder does not exist. Created.")  
    dir.create(foldername)
  }

  if("file_name" %ni% names(analysis_settings)){
    filename <- paste0(foldername,
                       analysis_settings[[2]],
                       "_",
                       analysis_settings[[3]],
                       "_",
                       analysis_settings[[4]],
                       "_",
                       analysis_settings[[5]],
                       "_",
                       analysis_settings[[6]],
                       extention
                      )
  }else if("file_name" %in% names(analysis_settings)){
    filename <- paste0(
                        foldername,
                        analysis_settings[[2]],
                        "_",
                        analysis_settings[[3]],
                        "_",
                        analysis_settings[[4]],
                        extention
                      )  
  }else{
    stop("Check the arguments of `analysis_settings`")  
  }

  return(filename)
}

save_figure <- function(p, analysis_settings, save_settings, filename=NULL){

  filename <- if_else(is.null(filename), make_filename(analysis_settings), filename)
  width <- save_settings$width
  height <- save_settings$height
  legend_ <- ifelse(is.null(save_settings$legend), TRUE, FALSE)
  themebw <- ifelse(is.null(save_settings$themebw), TRUE, FALSE)


  if(themebw)
    p <- p + theme_bw() 

  if(!legend_){
    p <- p + theme(legend.position="none", text = element_text(size=14), plot.title = element_text(hjust = 0.5))
  }else{
    p <- p + theme(text = element_text(size=14), plot.title = element_text(hjust = 0.5))
  }

  ggsave(filename, p, width=width, height=height)

  return(p)

}


check_figure_existence <- function(filename){
  # Check figure existence

  if(file.exists(paste0(filename))){
    return(TRUE)  
  }
  
  return(FALSE)

}


freq_parametric_figure <- function(res, analysis_settings, figure_settings){

  title <- if_else(is.null(figure_settings[["title"]]), "", figure_settings[["title"]])
  dodge <- if_else(is.null(figure_settings[["dodge"]]), -1/2, figure_settings[["dodge"]])
  
  if("Group" %in% colnames(res)){
    ggplot(res, aes(x=Model, y=estimate, colour=Group)) +
      geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
      geom_point(position = position_dodge(width = -1/2)) +
      geom_errorbar(size=0.5, width=0.3,
                    aes(ymin = estimate - 1.96 * std.error,
                                  ymax = estimate + 1.96 * std.error),
                    position = position_dodge(width = -1/2)
                    ) +
      coord_flip() +
      ggtitle(title) +
      xlab("Model") +
      ylab("Estimate") -> p

  }else{
    ggplot(res %>% filter(term!="(Intercept)"), aes(x=term, y=estimate, colour=Model)) +
      geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
      geom_point(position = position_dodge(width = -1/2)) +
      geom_errorbar(size=0.5, width=0.3,
                    aes(ymin = estimate - 1.96 * std.error,
                                  ymax = estimate + 1.96 * std.error),
                    position = position_dodge(width = -1/2)
                    ) +
      coord_flip() +
      ggtitle(title) +
      xlab("Covariate") +
      ylab("Estimate") -> p
  }

  return(p)

}


freq_nonparametric_figure <- function(res, analysis_settings, figure_settings){

  title <- if_else(is.null(figure_settings[["title"]]), "", figure_settings[["title"]])
  dodge <- if_else(is.null(figure_settings[["dodge"]]), -1/2, figure_settings[["dodge"]])
  
  if("Group" %in% colnames(res)){
    ggplot(res %>% filter(term %in% figure_settings$choices_rename), 
           aes(x=term, y=estimate, colour=Group)) -> p
  }else{
    ggplot(res %>% filter(term %in% figure_settings$choices_rename), 
           aes(x=term, y=estimate)) -> p
  }

  choices_rename <- figure_settings[["choices_rename"]]

  p +
    geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
    geom_point(position = position_dodge(width = -1/2)) +
    geom_errorbar(size=0.5, width=0.3,
                  aes(ymin = estimate - 1.96 * std.error,
                                ymax = estimate + 1.96 * std.error),
                  position = position_dodge(width = -1/2)
                  ) +
    scale_x_discrete(limits = rev(choices_rename)) +
    coord_flip() +
    ggtitle(title) +
    xlab("Covariate") +
    ylab("Estimate") -> p

  return(p)

}


stan_figure <- function(resStan, analysis_settings, figure_settings){

  rename_choices <- function(item, choices){
    n <- length(choices)

    if(item <= n)
      return(choices[as.integer(item)])
  }

  title <- if_else(is.null(figure_settings[["title"]]), "", figure_settings[["title"]])
  savefilename <- make_filename(analysis_settings)

  choices_rename <- figure_settings[["choices_rename"]]
  dodge <- if_else(is.null(figure_settings[["dodge"]]), -1/2, figure_settings[["dodge"]])

  focus_var_name <- figure_settings[["focus_var_name"]]
  focus_names <- analysis_settings[["interaction_name"]]

  if(is.null(focus_names)){
    result <- rstan::extract(resStan)
    values <- as_tibble((result$diff))
    values_qua <- data.frame(t(apply(values, 2, quantile, prob=c(0.05, 0.5, 0.95))))
    colnames(values_qua) <- c('p5', 'p50', 'p95')
    values_qua <- transform(values_qua, X=rownames(values_qua))

    values_qua %>%
      mutate(X = paste0("X", X)) %>%
      separate(X, into=c("temp", "Choice_int"), "V") %>%
      mutate(Choice_int = as.numeric(Choice_int)) %>%
      rowwise() %>%
      mutate(Choice = rename_choices(Choice_int, choices_rename)) %>%
      select(-Choice_int) -> values_qua

    values_qua$Choice <- factor(values_qua$Choice, level=choices_rename)

    p <- ggplot() +
              coord_flip() +
              geom_pointrange(data=values_qua, 
                              aes(x=Choice, y=p50, ymin=p5, ymax=p95), 
                              size=0.6) +
              scale_x_discrete(limits = rev(choices_rename)) +
              scale_y_continuous(labels=scales::percent_format()) +
              geom_hline(aes(yintercept=0), linetype="dashed") +
              xlab("Choice") + ylab("Estimate") + ggtitle(title)
  
  }else{
  
    focus_var_count <- length(focus_names)

    result <- rstan::extract(resStan)
    values <- as.tibble((result$diff))
    values_qua <- data.frame(t(apply(values, 2, quantile, prob=c(0.05, 0.5, 0.95))))
    colnames(values_qua) <- c('p5', 'p50', 'p95')
    values_qua <- transform(values_qua, X=rownames(values_qua))

    values_qua %>%
      separate(X, into=c("Choice_int", "Focus_int"), "\\.") %>%
      mutate(Choice_int = as.numeric(Choice_int), Focus_int =as.numeric(Focus_int)) %>%
      rowwise() %>%
      mutate(Choice = rename_choices(Choice_int, choices_rename),
             Group = rename_choices(Focus_int, focus_names)) %>%
      select(-Choice_int, -Focus_int) -> values_qua

    values_qua$Choice <- factor(values_qua$Choice, level=choices_rename)
    # values_qua$Group <- factor(values_qua$Group, level=focus_names)

    p <- ggplot() +
              coord_flip() +
              geom_pointrange(data=values_qua, 
                              aes(x=Choice, y=p50, ymin=p5, ymax=p95, colour=Group), 
                              size=0.6, position = position_dodge(width = -1/2)) +
              scale_x_discrete(limits = rev(choices_rename)) +
              scale_y_continuous(labels=scales::percent_format()) +
              geom_hline(aes(yintercept=0), linetype="dashed") +
              xlab("Choice") + ylab("Value") + ggtitle(title)
  
  }

  # Store values
  filename <- make_filename(analysis_settings, extention=".csv")
  values_qua %>%
    mutate(Table = paste0(
                          as.character(round(p50, 3)),
                          " [",
                          as.character(round(p5, 3)),
                          ", ",
                          as.character(round(p95, 3)),
                          "]"
                          )
          ) -> save_table
  write.csv(save_table, filename, row.names=F)

  return(p)

}


### Figure Finish}}}



### Analysis settings {{{

cov_types <- list(
                   # None = list(),
                   WithoutIncome = list(c("Q15.5", "age", "Age"), 
                                        c("Q15.2", "Q15.2_collapse", "Gender"),
                                        c("Q15.7", "Q15.7_collapse_new", "Education")
                                        )
                 )
cov_types_names <- names(cov_types)

interaction_types <- list(
                            None = list(),
                            IdealConst = list(c("Q14.3", "Near_collapse_new", "IdealConst")),  # (B) is 0
                            partyLDP = list(c("Party_LDP", "", "Party")),
                            partyIndp = list(c("Party_Indp", "", "Party")),
                            partyLeft = list(c("Party_Left", "", "Party")),
                            Knowledge = list(c("Q3_collapse", "", "Knowledge"))
                         )

interaction_names <- list(
                            None = list(),  # placeholder
                            IdealConst = list(c("Protect Rights", "Embody Ideals")),  # 0 to 1
                            partyLDP = list(c("Non-LDP", "LDP")),
                            partyIndp = list(c("Non-Independent", "Independent")),
                            partyLeft = list(c("Non-Left", "Left")),
                            Knowledge = list(c("Unknowledgeable", "Knowledgeable"))
                         )

analysis_functions <- list(
                            BayesSimple = analysis_bayes_ordered,
                            BayesSimpleFocus = analysis_bayes_ordered_focus,
                            BayesCov = analysis_bayes_cov,
                            BayesCovFocus = analysis_bayes_cov_focus,
                            FreqParametric = analysis_freq_parametric,
                            FreqNonparametric = analysis_freq_nonparametric
                          )


figure_functions <- list(
                            BayesSimple = stan_figure,
                            BayesSimpleFocus = stan_figure,
                            BayesCov = stan_figure,
                            BayesCovFocus = stan_figure,
                            FreqParametric = freq_parametric_figure,
                            FreqNonparametric = freq_nonparametric_figure
                          )


### }}}


### Final tile plot {{{

figure_all_choices <- function(analysis, selected = 1:length(use_order)) {
  
  choice <- c(1, 2, 3, 4, 5)
  
  map_dfr(selected,
          function(x) {
            right <- use_order[x]
            path <- paste0(savefolder_base, right, "/", sample_name, "_", right, analysis, ".obj")
            tmp <- readRDS(path)
            ggs_obj <- ggs(tmp$returned) %>% filter(grepl("diff", Parameter))
            geweke <- ggs_geweke(ggs_obj)$data
            geweke %>%
              filter(z < -2 | z > 2) %>%
              pull(Chain) %>% unique() -> exclude_chain
            
            tmp$returned %>%
              tidybayes::spread_draws(diff[choice]) %>%
              filter(! .chain %in% exclude_chain) %>%
              tidybayes::mean_hdi(.width = hdi_interval) %>%
              ungroup() %>%
              rename(point = diff) %>%
              mutate(type = "diff",
                     Right = right,
                     RightName = use_order_name[x],
                     choice = case_when(
                       choice == 1 ~ "Good",
                       choice == 2 ~ "Lean Good",
                       choice == 3 ~ "Neutral",
                       choice == 4 ~ "Lean Bad",
                       choice == 5 ~ "Bad"
                     )) -> a
            
            return(a)
          }
  ) -> res 
  
  res %>%
    mutate(choice = factor(choice, levels = c("Good", "Lean Good",
                                              "Neutral", "Lean Bad", "Bad")))  %>% 
    mutate(RightName = factor(RightName, levels = use_order_name)) -> res
 
  res %>%
    mutate(Category = case_when(grepl("Art 9", RightName) ~ "Article 9",
                                grepl("Emergency", RightName) ~ "National Emergency",
                                grepl("Family|Separation|voting", RightName) ~ "Cultural / Traditional",
                                grepl("Environment|Privacy|transparency", RightName) ~ "Progressive / Rights",
                                TRUE ~ "Else")) %>%
    mutate(Category = factor(Category, levels = c("Article 9", "National Emergency",
                                                  "Cultural / Traditional", "Progressive / Rights"))) -> res_use
  
  res_use %>%
    mutate(PointVal = scales::percent(point, accuracy = 0.1, suffix = "")) -> res_use
  
  res_use %>%  
    mutate(Significant = case_when(.lower < 0 & .upper < 0 ~ "A",
                                   .lower > 0 & .upper > 0 ~ "A",
                                   TRUE ~ "B")) -> res_use
  
  ggplot(res_use, aes(x = choice, y = point,
                      #linetype = Significant, colour = Significant, 
                      shape = choice)) +
    geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
    geom_linerange(aes(ymin = .lower, ymax = .upper), size = 0.6, position = position_dodge2(width = 1/3)) +
    geom_point(size = 2.2, position = position_dodge2(width = 1/3), fill = "white") +
    theme_bw() +
    facet_wrap(~ RightName, ncol = 3) +
    guides(linetype = FALSE,
           shape = FALSE,
           colour = FALSE) +
    scale_colour_manual(values = c("#1380A1", "gray40")) +
    scale_shape_manual(values = c(1, 19, 17, 15, 22)) +
    scale_y_continuous(limit = c(-0.1, 0.1), 
                       breaks = c(-0.1, -0.05, 0.0, 0.05, 0.1), labels = scales::percent_format(accuracy = 1)) +
    scale_linetype_manual(values = c("solid", "solid")) +
    #geom_text(aes(label = PointVal), position = position_nudge(x = 0.26)) +
    guides(shape = guide_legend(title = "Choice",
                                label.position = "top",
                                title.position = "left",
    )) +
    theme(legend.position = "top",
          axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.ticks.y = element_blank(),
          text = element_text(size = 13),
          strip.background = element_blank(),
          panel.border = element_rect(colour = "black"),
          panel.grid = element_blank()
    ) +
    xlab("") +
    ylab("Estimated effect") -> p ; p
  
  savepath <- paste0(savefolder_base, "result_tile/", "fig_all_choice", analysis, ".pdf")
  ggsave(filename = savepath, plot = p, width = 7, height = 8)
  
}


read_res <- function(analysis, selected = 1:length(use_order), factorize = TRUE) {
  
  choice <- c(1,2)
  
  map_dfr(selected,
          function(x) {
            right <- use_order[x]
            path <- paste0(savefolder_base, right, "/", sample_name, "_", right, analysis, ".obj")
            tmp <- readRDS(path)
            ggs_obj <- ggs(tmp$returned) %>% filter(grepl("diff", Parameter))
            geweke <- ggs_geweke(ggs_obj)$data
            geweke %>%
              filter(z < -2 | z > 2) %>%
              pull(Chain) %>% unique() -> exclude_chain
            
            tmp$returned %>%
              tidybayes::spread_draws(diff_D[choice]) %>%
              filter(! .chain %in% exclude_chain) %>%
              tidybayes::mean_hdi(.width = hdi_interval) %>%
              ungroup() %>%
              rename(point = diff_D) %>%
              mutate(type = "diff_D",
                     Right = right,
                     RightName = use_order_name[x]) -> a
            
            tmp$returned %>%
              tidybayes::spread_draws(diff_DwN[choice]) %>%
              #tidybayes::median_hdi(.width = hdi_interval) %>%
              tidybayes::mean_hdi(.width = hdi_interval) %>%
              ungroup() %>%
              rename(point = diff_DwN) %>%
              mutate(type = "diff_DwN",
                     Right = right,
                     RightName = use_order_name[x]) -> b
            
            return(bind_rows(a, b))
          }
  ) -> res 
  
  res %>%
    mutate(Choice = if_else(choice == 1, "Good", "Bad")) -> res
  
  if (!factorize) {
    return(res)
  }
  
  res %>%
    mutate(Choice = factor(Choice, levels = c("Bad", "Good"))) %>%
    mutate(RightName = factor(RightName, levels = use_order_name)) -> res
  return(res)
}


plot_res <- function(analysis, selected = 1:length(use_order)) {
  res <- read_res(analysis, selected)
  
  # p <- plot_ggplot(res, type = "diff_D") ; p
  # savepath <- paste0(savefolder_base, "result_tile/", "fig", analysis, "_diffD", ".pdf")
  # ggsave(filename = savepath, plot = p, width = 8, height = 5.5)
  #
  # p <- plot_ggplot(res, type = "diff_DwN") ; p
  # savepath <- paste0(savefolder_base, "result_tile/", "fig", analysis, "diffDwN", ".pdf")
  # ggsave(filename = savepath, plot = p, width = 8, height = 5.5)
  
  p <- plot_ggplot_2(res, type = "diff_D") ; p
  savepath <- paste0(savefolder_base, "result_tile/", "fig_single", analysis, "_diff", ".pdf")
  ggsave(filename = savepath, plot = p, width = 7, height = 6)
  
  # p <- plot_ggplot_2(res, type = "diff_DwN") ; p
  # savepath <- paste0(savefolder_base, "result_tile/GoodOnly/", "fig_single", analysis, "diffDwN", ".pdf")
  # ggsave(filename = savepath, plot = p, width = 7, height = 6)
}


plot_ggplot_2 <- function(res, type = c("diff_D", "diff_DwN")) {
  type <- match.arg(type)
  res_use <- res %>% filter(type == {{type}}) %>% 
              filter(Choice == "Good") %>%
              mutate(Significant = case_when(.lower < 0 & .upper < 0 ~ "A",
                                         .lower > 0 & .upper > 0 ~ "A",
                                         TRUE ~ "B"))
  
  res_use %>%
    mutate(Category = case_when(grepl("Art 9", RightName) ~ "Article 9",
                                grepl("Emergency", RightName) ~ "National Emergency",
                                grepl("Family|Separation|voting", RightName) ~ "Cultural / Traditional",
                                grepl("Environment|Privacy|transparency", RightName) ~ "Progressive / Rights",
                                TRUE ~ "Else")) %>%
    mutate(Category = factor(Category, levels = c("Article 9", "National Emergency",
                                                  "Cultural / Traditional", "Progressive / Rights"))) -> res_use
  
  res_use %>%
    mutate(PointVal = scales::percent(point, accuracy = 0.1, suffix = "")) -> res_use
  
  ggplot(res_use, aes(x = reorder(RightName, desc(RightName)), y = point,
                      linetype = Significant, colour = Significant, shape = Choice)) +
    geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
    geom_linerange(aes(ymin = .lower, ymax = .upper), size = 0.6, position = position_dodge2(width = 1/3)) +
    geom_point(size = 2.2, position = position_dodge2(width = 1/3), fill = "white") +
    coord_flip() +
    theme_bw() +
    facet_grid(Category ~ ., scales = "free_y") +
    guides(linetype = FALSE,
           shape = FALSE,
           colour = FALSE) +
    scale_colour_manual(values = c("#1380A1", "gray50")) +
    #scale_shape_manual(values = c(22, 19, 17, 15, 1)) +
    scale_y_continuous(breaks = c(-0.02, 0, 0.02, 0.04, 0.06, 0.08, 0.1), labels = scales::percent_format(accuracy = 1)) +
    scale_linetype_manual(values = c("solid", "solid")) +
    geom_text(aes(label = PointVal), position = position_nudge(x = 0.26)) +
    theme(legend.position = "top",
          axis.title.y = element_blank(),
          axis.ticks.y = element_blank(),
          text = element_text(size = 13),
          strip.background = element_blank(),
          panel.border = element_rect(colour = "black"),
          panel.grid = element_blank()
    ) +
    xlab("") +
    ylab("Estimated effect") -> p ; p
  return(p)
}


plot_ggplot <- function(res, type = c("diff_D", "diff_DwN")) {
  type <- match.arg(type)
  res_use <- res %>% filter(type == {{type}})
  
  text_df <- data.frame(Choice = c("Good", "Bad"),
                        point = c(-0.065, -0.065),
                        RightName = factor(c(use_order_name[1], use_order_name[1]), levels = use_order_name),
                        label = c("Good", "Bad"))
  
  ggplot(res_use, aes(x = Choice, y = point,
                      linetype = Choice, colour = Choice, shape = Choice)) +
    geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
    geom_linerange(aes(ymin = .lower, ymax = .upper), size = 0.6, position = position_dodge2(width = 1/3)) +
    geom_point(size = 2.2, position = position_dodge2(width = 1/3), fill = "white") +
    facet_wrap(~ RightName, ncol = 3) +
    coord_flip() +
    theme_bw() +
    geom_text(data = text_df, label = c("Good", "Bad")) +
    guides(linetype = FALSE,
           shape = FALSE,
           colour = FALSE) +
    scale_colour_manual(values = c("#990000", "#1380A1")) +
    scale_shape_manual(values = c(22, 19, 17, 15, 1)) +
    scale_y_continuous(labels = scales::percent) +
    scale_linetype_manual(values = c("twodash", "solid")) +
    theme(legend.position = "top",
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          text = element_text(size = 13),
          panel.grid = element_blank()
    ) +
    xlab("") +
    ylab("Estimated effect") -> p
}



read_res_focus <- function(analysis, selected = 1:length(use_order), focus_name = NULL, factorize = TRUE) {

  choice <- c(1,2)
  focus <- c(1,2)
  
  map_dfr(selected,
          function(x) {
            right <- use_order[x]
            path <- paste0(savefolder_base, right, "/", sample_name, "_", right, analysis, ".obj")
            tmp <- readRDS(path)
            ggs_obj <- ggs(tmp$returned) %>% filter(grepl("diff", Parameter))
            geweke <- ggs_geweke(ggs_obj)$data %>% mutate(pval = 2*(pnorm(-abs(z))))
            geweke %>%
              filter(z < -2 | z > 2) %>%
              pull(Chain) %>% unique() -> exclude_chain
            if (length(exclude_chain) == 5){
              exclude_chain <- NULL
            }
            
            tmp$returned %>%
              tidybayes::spread_draws(diff_D[choice, focus]) %>%
              filter(! .chain %in% exclude_chain) %>%
              tidybayes::mean_hdi(.width = hdi_interval) %>%
              ungroup() %>%
              rename(point = diff_D) %>%
              mutate(type = "diff_D",
                     Right = right,
                     RightName = use_order_name[x]) -> a
            
            tmp$returned %>%
              tidybayes::spread_draws(diff_DwN[choice, focus]) %>%
              filter(! .chain %in% exclude_chain) %>%
              tidybayes::mean_hdi(.width = hdi_interval) %>%
              ungroup() %>%
              rename(point = diff_DwN) %>%
              mutate(type = "diff_DwN",
                     Right = right,
                     RightName = use_order_name[x]) -> b
            
            return(bind_rows(a, b))
          }
  ) -> res 
  
  res %>%
    mutate(Choice = if_else(choice == 1, "Good", "Bad")) %>%
    mutate(Group = if_else(focus == 1, focus_name[1], focus_name[2])) -> res
  
  if (!factorize) {
    return(res)
  }
  
  res %>%
    mutate(Choice = factor(Choice, levels = c("Good", "Bad"))) %>%
    mutate(RightName = factor(RightName, levels = use_order_name)) %>%
    mutate(Group = factor(Group, levels = rev(focus_name))) -> res
  return(res)
}

plot_res_focus <- function(analysis, selected = 1:length(use_order), focus_name = NULL) {
  res <- read_res_focus(analysis, selected, focus_name)
  
  # p <- plot_ggplot_focus(res, type = "diff_D") ; p
  # savepath <- paste0(savefolder_base, "result_tile/", "fig", analysis, "_diffD", ".pdf")
  # ggsave(filename = savepath, plot = p, width = 8, height = 9)
  #
  # p <- plot_ggplot_focus(res, type = "diff_DwN") ; p
  # savepath <- paste0(savefolder_base, "result_tile/", "fig", analysis, "diffDwN", ".pdf")
  # ggsave(filename = savepath, plot = p, width = 8, height = 9)
  
  p <- plot_ggplot_focus_2(res, type = "diff_D") ; p
  savepath <- paste0(savefolder_base, "result_tile/", "fig_single", analysis, "_diff", ".pdf")
  ggsave(filename = savepath, plot = p, width = 11, height = 8)
  
  # p <- plot_ggplot_focus_2(res, type = "diff_DwN") ; p
  # savepath <- paste0(savefolder_base, "result_tile/GoodOnly/", "fig_single", analysis, "diffDwN", ".pdf")
  # ggsave(filename = savepath, plot = p, width = 11, height = 8)
}



plot_ggplot_focus_2 <- function(res, type = c("diff_D", "diff_DwN")) {
  type <- match.arg(type)
  res_use <- res %>% filter(type == {{type}}) %>%
              filter(Choice == "Good") %>%
                mutate(Significant = case_when(.lower < 0 & .upper < 0 ~ "A",
                                               .lower > 0 & .upper > 0 ~ "A",
                                               TRUE ~ "B"))
  
  res_use %>%
    mutate(Category = case_when(grepl("Art 9", RightName) ~ "Article 9",
                                grepl("Emergency", RightName) ~ "National Emergency",
                                grepl("Family|Separation|voting", RightName) ~ "Cultural / Traditional",
                                grepl("Environment|Privacy|transparency", RightName) ~ "Progressive / Rights",
                                TRUE ~ "Else")) %>%
    mutate(Category = factor(Category, levels = c("Article 9", "National Emergency",
                                                  "Cultural / Traditional", "Progressive / Rights"))) -> res_use
  
  res_use %>%
    mutate(PointVal = scales::percent(point, accuracy = 0.1, suffix = "")) -> res_use
  
  val_min <- min(res_use$.lower)
  val_max <- max(res_use$.upper)
  val_break <- seq(ifelse(val_min < -0.05, -0.1, -0.05), 
                   ifelse(val_max > 0.1, 0.15, 0.1), 
                   0.05)
  
  ggplot(res_use, aes(x = reorder(RightName, desc(RightName)), y = point,
                      colour = Significant)) +
    geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
    geom_linerange(aes(ymin = .lower, ymax = .upper, linetype = Significant), size = 0.6, position = position_dodge2(width = 1/3)) +
    geom_point(size = 2.2, position = position_dodge2(width = 1/3), fill = "white") +
    facet_grid(Category~ Group, scale = "free_y") +
    coord_flip() +
    theme_bw() +
    guides(shape = FALSE,
          linetype = FALSE,
          colour = FALSE) +
    scale_colour_manual(values = c("#1380A1","gray50")) +
    scale_shape_manual(values = c(19, 22, 17, 15, 1)) +
    scale_y_continuous(breaks = val_break, labels = scales::percent_format(accuracy = 1)) +
    scale_linetype_manual(values = c("solid", "solid")) +
    geom_text(aes(label = PointVal), position = position_nudge(x = 0.22)) +
    theme(legend.position="top",
          text = element_text(size = 13),
          strip.background = element_blank(),
          panel.border = element_rect(colour = "black"),
          panel.grid = element_blank()
    ) +
    xlab("") +
    ylab("Estimated effect") -> p ; p
}


plot_ggplot_focus <- function(res, type = c("diff_D", "diff_DwN")) {
  type <- match.arg(type)
  res_use <- res %>% filter(type == {{type}})
  
  val_max <- res_use %>% filter(Right == "A9DefenceAmy") %>% pull(.upper) %>% max() + 0.030
  text_df <- data.frame(Choice = c("Good", "Bad"),
                        point = c(val_max, val_max),
                        RightName = factor(c(use_order_name[1], use_order_name[1]), levels = use_order_name),
                        Group = c("", ""),
                        label = c("Good", "Bad"))
  
  ggplot(res_use, aes(x = Choice, y = point,
                      colour = Choice)) +
    geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
    geom_linerange(aes(ymin = .lower, ymax = .upper, linetype = Group), size = 0.6, position = position_dodge2(width = 1/3)) +
    geom_point(aes(shape = Group), size = 2.2, position = position_dodge2(width = 1/3), fill = "white") +
    facet_wrap(~ RightName, ncol = 3) +
    theme_bw() +
    geom_text(data = text_df, label = c("Good", "Bad")) +
    guides(shape = guide_legend(title = "Group",
                                label.position = "top",
                                title.position = "left"
    ),
    linetype = guide_legend(title = "Group",
                            label.position = "top",
                            title.position = "left"),
    colour = FALSE) +
    scale_colour_manual(values = c("#1380A1","#990000")) +
    scale_shape_manual(values = c(19, 22, 17, 15, 1)) +
    scale_y_continuous(labels = scales::percent) +
    scale_linetype_manual(values = c("solid", "twodash")) +
    theme(legend.position="top",
          axis.title.x=element_blank(),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          text = element_text(size = 13),
          panel.grid = element_blank()
    ) +
    xlab("") +
    ylab("Estimated effect") -> p
}

table_res <- function(type = c("diff_D", "diff_DwN"), significant_only = TRUE) {
  
  type <- match.arg(type)
  res <- read_res(analysis = "_WithoutIncome_None_BayesCov", factorize = FALSE)
  
  colname <- "ATE"
  res %>%
    filter(type == {{type}} & Choice == "Good") %>%
    { if (significant_only) filter(., (.lower < 0 & .upper <0) | (.lower > 0 & .upper > 0))
      else .} %>%
    select(RightName, point) %>%
    mutate(point = round(point * 100, digits = 1)) %>%
    mutate(type = colname) -> res_main
  
  res_list <- list()
  analyses <- c("_WithoutIncome_IdealConst_BayesCovFocus",
                "_WithoutIncome_partyLDP_BayesCovFocus",
                "_WithoutIncome_partyIndp_BayesCovFocus",
                "_WithoutIncome_partyLeft_BayesCovFocus")
  focus_names <- list(
    c("Protect Rights", "Embody Ideals"),
    c("Non-LDP", "LDP"),
    c("Non-Independent", "Independent"),
    c("Non-Left", "Left")
  )
  
  res_list[[1]] <- res_main
  j <- 2
  for (i in seq_along(analyses)) {
    res <- read_res_focus(analyses[i], focus_name = focus_names[[i]], factorize = FALSE)
    
    col_names <- unique(res$Group)
    res %>%
      filter(type == {{type}} & Choice == "Good") %>%
      { if (significant_only) filter(., (.lower < 0 & .upper <0) | (.lower > 0 & .upper > 0))
        else .} %>%
      select(RightName, point, Group) %>%
      mutate(point = round(point * 100, digits = 1)) %>%
      rename(type = Group) -> temp
    res_list[[j]] <- temp
    j <- j + 1
    
  }
  
  res_all <- bind_rows(res_list)
  res_all %>%
    {if (significant_only) add_row(., RightName = "Privacy", point = NA, type = "LDP") 
      else .} %>%
    spread(key = type, value = point) %>%
    left_join(x = tibble(RightName = use_order_name), y = .) %>%
    select(RightName, ATE, `Embody Ideals`, `Protect Rights`, LDP, `Left`, `Independent`) -> res_table
  
  if (significant_only) {
    savepath <- paste0(savefolder_base, "result_tile/", "table-", type, ".csv")
    write_csv(res_table, savepath) 
    
    p <- res_table_figure(res_table)
    savepath <- paste0(savefolder_base, "result_tile/", "table-fig-", type, ".pdf")
    ggsave(filename = savepath, plot = p, width = 10, height = 7.5)  
  } else {
    # show all values
    savepath <- paste0(savefolder_base, "result_tile/", "table-", type, "-all", ".csv")
    write_csv(res_table, savepath) 
    
    p <- res_table_figure(res_table)
    savepath <- paste0(savefolder_base, "result_tile/", "table-fig-", type, "-all", ".pdf")
    ggsave(filename = savepath, plot = p, width = 10, height = 7.5)
  }
}

res_table_figure <- function(res_table) {
  # Make it longer again
  res_table %>%
    gather(key = Type, value = point, -RightName) %>%
    mutate(Group = case_when(Type == "ATE" ~ "No interaction",
                             TRUE ~ "Interaction" )) %>% 
    mutate(Group = factor(Group, levels = c("No interaction", "Interaction")),
           RightName = factor(RightName, levels = use_order_name)) %>%
    mutate(Type = factor(Type, levels = c("ATE", "Embody Ideals", "Protect Rights", "LDP", "Left", "Independent"))) %>%
    mutate(Category = case_when(grepl("Art 9", RightName) ~ "Article 9",
                                grepl("Emergency", RightName) ~ "National Emergency",
                                grepl("Family|Separation|voting", RightName) ~ "Cultural / Traditional",
                                grepl("Environment|Privacy|transparency", RightName) ~ "Progressive / Rights",
                                TRUE ~ "Else")) %>%
    mutate(Category = factor(Category, levels = c("Article 9", "National Emergency",
                                                  "Cultural / Traditional", "Progressive / Rights"))) -> res_heat
  
  ggplot(res_heat, aes(x = Type, y = reorder(RightName, desc(RightName)), fill = point)) +
    geom_tile() +
    scale_fill_gradient(name = "Effect size", na.value = 'white', low = "gray90", high = "blue") +
    facet_grid(Category ~ Group, scale = "free", space = "free") +
    xlab("Type of effect") +
    ylab("") +
    theme_bw() +
    theme(legend.position="top",
          text = element_text(size = 13),
          strip.background = element_blank(),
          panel.border = element_rect(colour = "black"),
          panel.grid = element_blank()
    ) -> p ; p
}

### }}}


## Descriptive stats
data_description <- function(data) {
  # Knowledge score
  data$Q3.2 <- apply(data[, "Q3.2"], 1,
                     function(x){
                       change_choices(if_else(!is.na(x), x, "NA"),
                                      "Q3.2") 
                     })
  
  data$Q3.3 <- apply(data[, "Q3.3"], 1,
                     function(x){
                       change_choices(if_else(!is.na(x), x, "NA"),
                                      "Q3.3") 
                     })
  
  data$Q3.4 <- apply(data[, "Q3.4"], 1,
                     function(x){
                       change_choices(if_else(!is.na(x), x, "NA"),
                                      "Q3.4") 
                     })
  
  data %>%
    rowwise() %>%
    mutate(Q3_CorrectNum = sum(c(Q3.2=="Correct", Q3.3=="Correct", Q3.4=="Correct"))) %>%
    mutate(Q3 = case_when(
      Q3_CorrectNum == 3 ~ 1,   
      Q3_CorrectNum == 2 ~ 1,   
      TRUE ~ 0  
    )
    ) -> data 
  
  # Party ID
  sapply(1:nrow(data), function(i) {
    if (data[i, "Q13.2"] == "どの政党も支持しない") {
      return(data[i, "Q13.3"] %>% as.character())
    } else {
      return(data[i, "Q13.2"] %>% as.character())
    }
  }) -> data$PartyID
  
  
  data$Party_LDP <- apply(data[, c("Q13.2", "Q13.3")], 1,
                          function(x){
                            change_choices(x, "party_LDP_collapse")  
                          })
  
  data$Party_Indp <- apply(data[, c("Q13.2", "Q13.3")], 1,
                           function(x){
                             change_choices(x, "party_indp_collapse")  
                           })
  
  data$Party_Left <- apply(data[, c("Q13.2", "Q13.3")], 1,
                           function(x){
                             change_choices(x, "party_left_collapse")  
                           })
  
  
  # Ideal
  data$Q14.3 <- apply(data[, "Q14.3"], 1,
                     function(x){
                       change_choices(if_else(!is.na(x), x, "NA"),
                                      "Near_collapse_new") -> val
                       if (val == 0) 
                         return("Protect Rights")
                       if (val == 1)
                         return("Embody Ideals")
                       return("NA")
                     })
  
  return(data)
} 


check_rhat <- function() {
  mcmc_files <- list.files(path = "figures300s/", pattern = "BayesCov.*\\.obj",
                           full.names = TRUE, recursive = TRUE)
  
  for (path in mcmc_files) {
    stan_obj <- readRDS(path)$returned
    ggs_obj <- ggs(stan_obj) %>% filter(grepl("diff", Parameter))
    Rhat <- ggs_Rhat(ggs_obj)$data$Rhat
    if (sum(Rhat > 1.1) != 0)
      print(path)
  }
}
